Compare NEXRAD Radar with KAZR#

Within this notebook, we will look at comparing NEXRAD data with data collected from the KAZR radar (vertically pointing)

Data#

We have two data sources:

NEXRAD Data#

NEXRAD data is stored on AWS through the Open Data Initiative

There is a specific Python package called nexradaws, which makes working with these data easier!

The nearest NEXRAD radars to the TRACER field campaign are:

We will be accessing data on the NEXRAD AWS Bucket (maintained by Unidata)

KAZR Data#

The KAZR data explored here is exploratory and has not been calibrated, indicating that it is not meant to be used for analysis, but rather a proof-of-concept. It is a vertically pointing Ka-band radar.

From the official ARM KAZR Documentation, “The Ka-band ARM zenith radar (KAZR) remotely probes the extent and composition of clouds at millimeter wavelengths. The KAZR is a zenith-pointing Doppler radar that operates at a frequency of approximately 35 GHz. The main purpose of this radar is to determine the first three Doppler moments (reflectivity, vertical velocity, and spectral width) at a range resolution of approximately 30 meters from near-ground to nearly 20 km in altitude.”

After finding data on the ARM Data Hub using this query, download the data and place in some data directory!

Imports#

import xarray as xr
import pyart
import nexradaws
import tempfile
from datetime import datetime
import pytz
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from distributed import LocalCluster, Client
import numpy as np
import xcollection as xc
import cartopy.feature as cfeature
import hvplot.xarray
import holoviews as hv
import warnings
warnings.filterwarnings("ignore")
hv.extension('bokeh')

cluster = LocalCluster()
client = Client(cluster)
client

2022-03-23 12:35:37,240 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-752c9wb7', purging
2022-03-23 12:35:37,242 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-fy_90cvl', purging
2022-03-23 12:35:37,242 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-wc05egad', purging
2022-03-23 12:35:37,243 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-yjd669ot', purging
2022-03-23 12:35:37,243 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-4odlboz9', purging
2022-03-23 12:35:37,244 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-ib4c6d9k', purging
2022-03-23 12:35:37,244 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-_c3t0573', purging
2022-03-23 12:35:37,245 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-buisdlqx', purging
2022-03-23 12:35:37,245 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-gt_1tup4', purging
2022-03-23 12:35:37,246 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-_9y7swk5', purging
2022-03-23 12:35:37,246 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-2ixfjb7f', purging
2022-03-23 12:35:37,247 - distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/tracer-radar-analysis/notebooks/march-22-2022/dask-worker-space/worker-6a7n_3s5', purging

Client

Client-a6e75608-aacf-11ec-95bb-acde48001122

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

Access NEXRAD Data#

We can use the NexradAwsInterface class from nexradaws to acces the data, using the following steps:

  • Establish an initial connection

  • Query/subset for our desired site/time

  • Download to some directory

#https://stackoverflow.com/questions/32237862/find-the-closest-date-to-a-given-date/32237949
def nearest(items, pivot):
    return min(items, key=lambda x: abs(x - pivot))

def download_radar_data(site='KLOT',
                        cache_location=tempfile.mkdtemp(),
                        target_time=datetime.utcnow()):
    """
    Accesses data from AWS S3
    
    Parameters
    ----------
    site: str, default='KLOT' (Chicago)
       Site to download data from (ex. 'KLOT')
    cache_location: str, default=tempfile.mkdtemp()
       Where to place the cached data files
    target_time: datetime.datetime, default=datetime.utcnow()
       Time index to search for
    
    Returns
    -------
    file : str
        Path to the NEXRAD file
    """
    
    # Convert the datetime to UTC
    target_time = target_time.astimezone(pytz.utc)
    
    # Setup the initial connection
    aws_connection = nexradaws.NexradAwsInterface()
    
    # Check for available scans
    scans = aws_connection.get_avail_scans(target_time.year,
                                           target_time.month,
                                           target_time.day,
                                           site)
    
    # Retrieve the scan times
    scan_times = [scan.scan_time for scan in scans]
    
    # Loop through and find good scans
    good_scans = []
    good_times = []
    for i in range(len(scans)):
        if scan_times[i] is not None:
            good_times.append(scan_times[i])
            good_scans.append(scans[i])

    print(f'Found {len(good_scans)} good scans')
    print(f'Found {len(good_times)} good times')
    
    # Find the nearest time using our helper function
    nearest_time = nearest(good_times, target_time)
    
    # Locate closest time in our scans
    closest_time_index = good_times.index(nearest_time)
    
    # Download the file
    file = aws_connection.download(good_scans[closest_time_index],
                                   f'{cache_location}/{site}')
    
    # Return the filepath
    return file.success[0].filepath

Download the file#

time = datetime(2022, 3, 21, 7, 55, 0)
KLCH_file = download_radar_data(site='KLCH',
                                cache_location='/Volumes/Samsung_T5/argonne/nexrad-data',
                                target_time=time)

KHGX_file = download_radar_data(site='KHGX',
                                cache_location='/Volumes/Samsung_T5/argonne/nexrad-data',
                                target_time=time)

Found 235 good scans
Found 235 good times
Downloaded KLCH20220321_125135_V06
1 out of 1 files downloaded...0 errors
Found 296 good scans
Found 296 good times
Downloaded KHGX20220321_125323_V06
1 out of 1 files downloaded...0 errors

Run Through the Analysis#

Now that we have the data, let’s start our analysis!

Convert to an Xarray Object#

We add a couple of helper functions to convert each sweep to a dataset, then the full volume to a dictionary of datasets

def create_dataset_from_sweep(radar, sweep=0, field=None):
    """
    Creates an xarray.Dataset from sweeps
    """
    
    # Grab the range
    range = radar.range['data']
    
    # Grab the elevation
    elevation = radar.get_elevation(sweep)
    
    # Grab the azimuth
    azimuth = radar.get_azimuth(sweep)
    
    # Grab the x, y, z values
    x, y, z = radar.get_gate_x_y_z(sweep)
    
    # Grab the lat, lon, and elevation
    lat, lon, alt = radar.get_gate_lat_lon_alt(sweep)
    
    # Add the fields
    field_dict = {}
    for field in list(radar.fields):
        field_dict[field]=(["azimuth", "range"], radar.get_field(sweep, field))
    
    ds = xr.Dataset(
        data_vars=field_dict,
         coords=dict(
             x=(["azimuth", "range"], x),
             y=(["azimuth", "range"], y),
             z=(["azimuth", "range"], z),
             lat=(["azimuth", "range"], lat),
             lon=(["azimuth", "range"], lon),
             alt=(["azimuth", "range"], alt),
             azimuth=azimuth,
             elevation=elevation,
             sweep=np.array([sweep]),
             range=range),
    )
    
    return ds.chunk()

Create an “Xradar” object#

For now, we utilize xcollection which is a dictionary of datasets with some added functionality

def convert_to_xradar(radar):
    """
    Converts from radar to xradar
    """
    ds_list = []
    sweeps = radar.sweep_number['data']
    for sweep in sweeps:
        ds_list.append(create_dataset_from_sweep(radar, sweep))
    
    # Convert the numpy array to a list of strings
    dict_keys = [x for x in list(sweeps.astype(str))]
    
    # Zip the keys and datasets and turn into a dictionary
    dict_of_dsets = dict(zip(dict_keys, ds_list))
    
    return xc.Collection(dict_of_dsets)

Run our workflow#

Read the data into PyART#

klch_radar = pyart.io.read(KLCH_file)
khgz_radar = pyart.io.read(KHGX_file)

klch_radar.time['units']

'seconds since 2022-03-21T12:51:35Z'

Convert to Xradar#

klch_xradar = convert_to_xradar(klch_radar)
khgz_xradar = convert_to_xradar(khgz_radar)

klch_xradar['0']

klch_xradar['0']

<xarray.Dataset>
Dimensions:                       (azimuth: 720, range: 1832, elevation: 720,
                                   sweep: 1)
Coordinates:
    x                             (azimuth, range) float64 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    y                             (azimuth, range) float64 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    z                             (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    lat                           (azimuth, range) float64 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    lon                           (azimuth, range) float64 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    alt                           (azimuth, range) float64 dask.array<chunksize=(720, 1832), meta=np.ndarray>
  * azimuth                       (azimuth) float64 7.251 7.748 ... 6.254 6.748
  * elevation                     (elevation) float32 0.5136 0.5466 ... 0.4834
  * sweep                         (sweep) int32 0
  * range                         (range) float32 2.125e+03 ... 4.599e+05
Data variables:
    differential_phase            (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    spectrum_width                (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    cross_correlation_ratio       (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    clutter_filter_power_removed  (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    differential_reflectivity     (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    velocity                      (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>
    reflectivity                  (azimuth, range) float32 dask.array<chunksize=(720, 1832), meta=np.ndarray>

Visualize our Output#

Start first with xarray! We use matplotlib directly, specifying our axes

fig = plt.figure(figsize=(10,8))
ax = plt.subplot(111, projection=ccrs.PlateCarree())

khgz_xradar['0'].reflectivity.plot(ax=ax,
                                   y='lat',
                                   x='lon',
                                   cmap='pyart_HomeyerRainbow',
                                   vmin=-20,
                                   vmax=60)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.STATES)

<cartopy.mpl.feature_artist.FeatureArtist at 0x195b95730>
../_images/compare-kazr-nexrad_23_1.png 
khgz_xradar['0'].reflectivity.hvplot.quadmesh(x='lon',
                                              y='lat',
                                              rasterize=True,
                                              width=600,
                                              height=400,
                                              clim=(-20, 60),
                                              projection=ccrs.PlateCarree(),
                                              cmap='pyart_HomeyerRainbow',
                                              features=['land', 'lakes','borders','ocean'])